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Summary 

A general approach to the construction of convergence acceleration methods for vector sequences 
is proposed. Using this approach, one can generate some known methods, such as the minimal 
polynomial extrapolation, the reduced rank extrapolation, and the topological epsilon algorithm, 
and also some new ones. Some of the new methods are easier to implement than the known methods 
and are observed to have similar numerical properties. The convergence analysis of these new 
methods is carried out, and it is shown that they are especially suitable for accelerating the 
convergence of vector sequences that are obtained when one solves linear systems of equations 
iteratively. A stability analysis is also given, and numerical examples are provided. The convergence 
and stability properties of the topological epsilon algorithm are likewise given. 


1. Introduction 

In a recent work (ref. 1) a survey of convergence acceleration methods for sequences of vectors is 
given, and five of these methods are tested and compared numerically using a process that has been 
termed cycling : the minimal polynomial extrapolation (MPE), the reduced rank extrapolation 
(RRE), the scalar epsilon algorithm (SEA), the vector epsilon algorithm (VEA), and the topological 
epsilon algorithm (TEA). One of the conclusions of this survey is that the MPE and the RRE have 
about the same properties, and in general, have better convergence than others, in the sense that the 
MPE and the RRE achieve a given level of accuracy with fewer vectors than the SEA, VEA, and 
TEA. VEA and SEA are also similar in performance, except that the latter is more prone to 
numerical instability problems. However, TEA, while interesting from a theoretical point of view, 
appears to be not as effective as either VEA or SEA; see reference 1 for further details. 

All of the methods above have the following important properties: 

(1) It is observed numerically that in many instances they accelerate the convergence of a slowly 
converging vector sequence and they make a diverging sequence converge to an “anti-limit” that has 
an immediate interpretation. 

(2) They depend solely on the given vector sequence whose convergence is being accelerated; they 
do not depend on how the vector sequence is generated. 

(3) Their implementation is straightforward. 

For more details and an extensive bibliography see reference 1 . 

It turns out that the implementation of the MPE and RRE requires the least -squares solution of an 
overdetermined and in general inconsistent set of linear equations, the number of the equations in 
this set being equal to the dimension of the vectors in the given sequence. For many practical 
problems, the dimension of these vectors may be finite but very large; consequently, one may have to 
store a large rectangular matrix in memory, making the MPE and RRE somewhat expensive in both 
storage and time. Therefore, it would be desirable to have methods as efficient as MPE and RRE but 
less demanding in storage and time. 

In the next section a general framework for deriving convergence acceleration methods for vector 
sequences - vector accelerators for short - is proposed. Within this framework one can derive several 
methods, some old (MPE, RRE, and TEA) and some new. It turns out that one of the new methods is 
very similar to the MPE and RRE, but does not require the use of the least-squares method, requires 
very little storage, and, at least numerically, is as efficient as MPE. The convergence analysis of this 
method, which we shall call the modified MPE, (MMPE), is carried out in Section 3 for a class of 
vector sequences that includes those arising from the iterative solution of systems of linear equations. 
We prove that this method is a bona fide convergence acceleration method, and its rate of 
acceleration is also provided. The stability properties of this method are taken up in Section 4. In 
Section 5 the convergence and stability properties of the TEA are analyzed using the techniques of 
Sections 3 and 4. Finally in Section 6 we test the MMPE on some examples numerically, and compare 
it with the MPE. On the basis of this comparison one could conclude that the MPE and MMPE have 
similar performances, although there is as yet no complete theory to support this empirical 
conclusion. 


2. Development of Vector Accelerators 

In this section we shall develop a general framework within which one can derive vector 
accelerators of different kinds. We shall motivate this development in the way Shanks (ref. 2) 
motivates his development of the ^-transformation for scalar sequences. 

2.1. The Shanks Transformation 

Shanks starts with a scalar sequence X m , m = 0, 1, . . ., that has the property 


X m ~S+ £ a/XJ" as m — oo, (2.1) 

i=i 

where S, a it and X,- are constants independent of m, X,^l, /= 1, 2, . . ., X,^Xy, for all i^j, and 
|Xi|s|X 2 |a:. ... In (2.1) S is lim X m if all |X,-|< 1; otherwise, S is called the anti-limit of the sequence 

/W — 00 

X m , m = 0, 1, .... As one way of approximating 5, Shanks proposes to solve the set of 2£+l 
nonlinear equations 

k 

x m = s n,k+ £ ofkf, n<m<n + 2k, (2.2) 

/=i 

for S n j c , which is taken to be an approximation to S, with a,, X/, t = l k, being the rest of the 

unknowns. The solution turns out to have the following determinant representation: 



X n+1 ... 

X n + k 

AX n 

A X n + 1 . . . 

&Xn+k 

AXn+k-l 

AX n +k • * * 

&Xn + 2k- 

1 

i 

1 

AX„ 

AX n + l . . . 

A X n+k 


&Xn + k-l &Xn + k ••• ^n + 2k-\ 


where A is the forward difference operator defined by Ab,=b ; + j — bj, A p bj = A(AP ~ p> 2, 
provided the determinant in the denominator of (2.3) is nonzero. 

Two equivalent formulations follow from (2.3), independent of any reference to the nonlinear 
equations in (2.2). 

(a) With * as given in (2.3), there are k parameters /3„ /= 0, 1, . . ., k- 1, which solve the 
system of k+ 1 linear equations 


n <m <n + k. 


(2.4) 


k — l 

X m — Sn, k + 2 PiAX m + /, 

/ — 0 

as can be verified by solving (2.4) for S n tk by Cramer’s rule. By taking the differences of the 
equations in (2.4), we see that the j 8 / satisfy 


k— 1 

AX m = 2 0iA 2 X m +b n<m<n + k- 1, (2.5) 

i=0 

and that, once the ft* have been determined from (2.5), can be computed from one of the 
equations in (2.4), say that for which m=n. 

(b) From (2.4) it follows that S n>fc , along with the parameters y i9 / = 0, 1, k, satisfies the system 
of linear equations 

k 

s n ,k= Y) yi x m+i> n<m<n + k, ( 2 . 6 ) 

/' = 0 

subject to 

in=i- w 

/ = 0 

By taking the differences of the equations in (2.6), we see that the 7 ; satisfy 

k 

0=2 yiAXm + i* n<m<n + k- 1 , ( 2 . 8 ) 

/=o 

subject to (2.7). Once the 7 / have been determined from (2.7) and (2.8), S ntk can be computed from 
one of the equations in (2.6). Furthermore, if 7 * 5 * 0 , then (2.7) and (2.8) are equivalent to 

k - 1 

2 CjAX m + j = - AX m+k , n<m<n + k- 1, (2.9) 

/=o 

where 


7 / = — , 0 </<£, c k — 1 , 

k 

Lcj 

j= 0 


( 2 . 10 ) 


k 

provided 2 9 5 * 0 . 

y-o 

It has been proved by Wynn (ref. 3) that S ntk , when computed from a sequence X m , m = 0, 1 , . . ., 
that is of the form given in (2. 1 ), converges to S as n — 00 (k fixed), under certain conditions on the X* 
faster than X n itself. Wynn actually gives rates of convergence for S ntk for 00 . 
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2.2. Derivation of Vector Accelerators 

Let us now consider a sequence of vectors, x m , m = 0, 1, . . ., in a general normed vector space B, 
satisfying 


00 

x m ~s+ £ v i^T asm ~°° (2.11) 

i=i 

where s and t>/ are vectors in B, and X,- are scalars, independent of m, X,^l,/=1,2, . . ., A, Ay for all 
iytj, and (A^IA^ .... We also assume that there can be only a finite number of A,- having the same 
modulus. The meaning of (2.11) is that, for any integer N> 0, there exist a positive constant K and a 
positive integer rn 0 that depend only on N, such that for every m > wq 


N— 1 

ii (2-Ha) 

/= i 

with || * || being the norm associated with the vector space B. 

A simple example of such a sequence is that produced by a matrix iterative technique for solving 
the linear system of equations 

x=Ax+b (2.12) 

where A is a nondefective MxM matrix, and b and x are M-dimensional column vectors. If 5 is the 
solution to (2.12), and for given x 0 , the vectors x m are generated by 

x m + i=Ax m + b, m = 0,1,..., (2.13) 


then 

AT 

•*m = s+ £ otpiX 1 ?, m= 0,1,... (2.14) 

»=l 

where are scalars, X/ and v x are the eigenvalues and corresponding eigenvectors of the matrix A, 
and M' < M is the number of the distinct eigenvalues. 

The condition stated in (2.11) is analogous to that stated in (2.1) for scalar sequences. Since the 
Shanks transformation accelerates the convergence of scalar sequences satisfying (2.1), we expect 
that its extensions to vector sequences, through the formulations (a) and (b) following (2.3), may also 
produce acceleration of convergence for vector sequences satisfying (2. 1 1). The extensions of the two 
formulations can be achieved as follows: 

Approach (a) 

In equations (2.5) replace Xjbyxj, and ‘’solve” in some sense the resulting overdetermined (and in 
general inconsistent) system 


= £ &A lx m+i. n<m<n + k~ 1, (2.15) 

»=o 

for the j8/. Once the j3/ have been determined, compute the approximation s n ^ to 5 by 


k - 1 

s n,k =x n ~ ^ 0iAX n +j, 
i—0 


(2.16) 
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which is obtained by replacing S n>k and Xj in (2.4) by s n>k and Xj, respectively, and considering m = n. 

Approach (b) 

In equations (2.9) replace Xj by Xj, and “solve” in some sense the resulting overdetermined (and in 
general inconsistent) system 


k-\ 

Y c,^x m + i = -Ax m+k , n<m<n + k- 1, (2.17) 

( = 0 

for the Cj. Once the c it / = 0, 1 k- 1, have been determined, set 

It c*=l, 7 /= — — — , 0<j<k, (2.18) 

k 

E c j 

j = o 

k 

assuming ^ cy s*0. Finally compute the approximation s nt k to s by 

7 — 0 
k 

s n,k~ ^ Ji x n + i> ( 2 * 19 ) 

/ = 0 

which is obtained by replacing S ntk and Xj in (2.6) by and xj, respectively, and considering m=n. 

We see that for both approaches, we need to “solve” an overdetermined and, in general, 
inconsistent system of equations of the form 

k - 1 

Y djw m+i = w m , n <m <n + k- 1, (2.20) 

i = 0 

where wy and u>y are vectors in B , and dj are unknown scalars. If r, the dimension of B , is greater than 
k t then even one of the equations in (2.20) gives rise to an overdetermined system of r equations. We 
can, however, propose various ways for obtaining a set of d\ that solves (2.20) in some sense. In what 
follows, we give three such methods, with the understanding that other methods can also be 
proposed. 

Method (1) 

Assuming r> k 3 consider only one of the equations in (2.20), namely that with m = n, and solve for 

k - 1 

the dj that minimize some norm of the vector A = ^ d{w n+i — w n . Depending on the norm being 

/— 0 

used, different acceleration methods will be obtained. For example, if r is finite and the (weighted) t p 
norms are used with p= 1, 2, oo, then the determination of the rf,- becomes relatively easy. For p = 2 
the solution can be achieved by using any one of the least-squares packages available, and for /?= 1 
and/?= oo the minimization problems can be solved by using linear programming techniques; see the 
^ review paper by Rabinowitz (ref. 4). The f 2 norm with equal weights gives rise to RRE for Approach 

(a) and MPE for Approach (b). The rest of the acceleration methods have not appeared in the 
literature before. 

Method (2) 

Assuming r>k t consider only one of the equations in (2.20), namely that with m = n 3 and obtain 
the df by solving the system of k equations 
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( 2 . 21 ) 


k — 1 

L diQj(w„ +i ) =Qj(w„), j= 1 , . . .,k, 

;= 0 

where Qj are linearly independent bounded linear functionals on the space When 5 is an inner 
product space, we can take Qj{y) = (Qj,y) where qj are vectors in B and (v) is the inner product 
associated with B . If r is finite, and the vector qj is chosen to be the ij th unit vector, that is, 
{qj, z) = ij ih component of z, then the method above is equivalent to demanding that only k out of the 
r equations be satisfied, namely those corresponding to the /) th components, y = 1, . . k. Obviously 
such an acceleration method demands less storage and time for its implementation than methods like 
the MPE and RRE. It is not difficult to see that Approaches (a) and (b) both give the same 
acceleration method, which has not been given in the literature before. Due to its similarity to the 
MPE, we shall call this method the modified MPE (MMPE). In Sections 3 and 4 we shall analyze the 
convergence and stability properties of this method in detail. 

Method (3) 

Consider all the equations in (2.20) and obtain the d} by solving the system of k equations 


k- 1 

£ diQ(w m+i ) =Q(w m ), n<m<n + k-\, (2.22) 

/= o 

where Q is a bounded linear functional on the space B . In this case Approaches (a) and (b) give the 
same method, and this method is nothing but the TEA. In Section 5 we shall analyze the convergence 
and stability properties of this method in detail. By comparing (2.21) and (2.22), we see that for k= 1 
the MMPE and the TEA are identical when we choose Q\ = Q . 

Finally we note that all of the methods obtained as above are nonlinear in the 


3. Convergence Analysis Of MMPE 

As in Section 2.2, we start with a sequence of vectors x if / = 0, 1, . . ., in a normed vector space B 
with norm || • || , that has a limit or anti-limit s. We write W/= Ax/=X/+i — / = 0, 1, . . . . Then the 

MMPE, as obtained from Approach (a) or (b) in conjunction with Method (2) (see Section 2.2), can 
be summarized (and reformulated) as follows: By Approach (b), the approximation to s is given 
as 


k 

s n,k~ ^ 7 i x n+i> ( 3 * 1 ) 

/ = 0 

where 7 / are obtained from 
k 

L 7/ = 1 
1=0 

(3.2) 

k 

L r , QMn +i ) =0, 1 < j < k . 

i = 0 

When 7 * 5 * 0 , equations (3.2) are equivalent to (2.18) and (2.21), in which di=c iy vv w + / =w„ +/ *, 
0 </<£-!, and w n = as can be verified by inspection. 
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We denote the scalars Qj(,u m ) by u m j for 1 <j<k and m> 0, and we define £>(ct q , ctj, . . o k ) to 
be the determinant 


°0 

<7i 

O k 

“n, 1 

M n+l,l • 

• • u n + k, 1 

2 

u n + 1,2 • 

• • + 2 


D(a 0 , aj, . . o k ) = 


u n,k u n + l,k • • • u n + k,k 


(3.3) 


when <jj are scalars. Let N, be the cofactor of u, in the first row expansion of this determinant. Then 
k 

D(a 0 , ffi, . . ., a k ) = OjNi. (3.4) 

i = 0 

When <7/ are vectors, we again let £>(<7 0 , a u . . ., a k ) be defined by the determinant in (3.3), and take 
(3.4) as the interpretation of this determinant. Thus £>(<7 0 , <7i, . . ., a k ) is a scalar (or vector) if the <7, 
are scalars (or vectors). 

Solving the system in (3.2) by Cramer’s rule, we obtain 


Ni Ni 

7/ = =~ k . 0</<A:, 

1 1) £ Nj 

7 = 0 

provided £>( 1, 1, . . 1) 5*0; in what follows, we assume that this is so. 

By (3.5), (3.1), and (3.4), we can finally express s nt k as 

_ D ( X n , X n + i , . • • > *n + k) 

S "- k ~ D( 1, 1, . . ., 1) 


Lemma 3.1: The error in s n can be expressed as 


Sn,k~ s = 


D {x n s f x n ^ | s, 
D( 1,1,. 


•> x n+k s) 

1 ) 


(3.5) 


(3.6) 


(3.7) 


Proof: (3.7) follows easily from (3.4) to (3.6). 

□ 

In the sequel we shall assume that the vectors x m satisfy (2.1 1). Without loss of generality we shall 
also assume that A/ 5*0 and t>/5*0 for all / >1. Then 


u m ~ Z/A™ asm — oo, (3.8) 

i=l 
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where z,= (X,- l)u,-, /'= 1, 2 Since X,^l and t;,^0 for all /si, we have z,^0 for all /si. In 

addition, by (2.11), for any operator T in the dual space of B, 


T(u m )~ £ r(z/)XJ" as m— oo. (3.9) 

1=1 

Consequently 


~ asm-* 00, (3.10) 

/=1 

where z/j= Gy (z/), i^l, l<y<A:. 

Note that when the sequence x m , m- 0, 1, . . is generated by the matrix iterative method 
described in Section 2.2, the summations over / in (3.8), (3.9), and (3.10) extend as far as AT, which is 
a finite number; therefore, (3.9) and hence (3.10) hold automatically for this case, and — is replaced 
by =. 

The following theorem is the main result of this section. 

Theorem 3.2: Define Qj{vj) =Vjj, / > 1 , 1 <y <fc, and let 


V\,\ 

”2,1 • 

• ”*,1 

v l,2 

”2,2 • 

• v k,2 

V\,k 

v 2,k • ■ 

• v k,k 


(3.11) 


Assume that the u/ are linearly independent, and that the X/ satisfy 

|Xi|s. . .s|X^>|X Ar + i|s|X A . +2 |>. . .. (3.12) 

Then, for all sufficiently large n, D( 1, 1, . . ., 1)5*0; hence, s nt k , as given in (3.6), exists. 
Furthermore, 

s»,*-'S=r(/i)X2 +1 [l+o(l)]as«-oo, (3.13) 

where the vector T(/i) is nonzero and bounded for all sufficiently large n. If, in addition, 
l\t+il>l\t + 2 l> then 


”l 

”2 • 

•• v k+ 1 




*1.1 

*2,1 • 

• • Zk+ 1,1 




*1,2 

*2,2 • 

• • ** + 1,2 

k 

n 

/ = i 

(\t+l - M 
(X,-l)2 • 

(3.14) 

Zl,k 

*2,k ■ ■ 

■ • Zk+\,k 
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Proof: For simplicity of notation we shall sometimes denote G n =D(x n —s, . . x n+k —s) and 

H n =D( 1 1), and we shall shorten “ot n ~l3 n as n — oo” to 

By (3.3) and (3.10) we have 


H n ~ 


oo 

OO 


oo 


£ 


• £ 

'1 = 1 

'1 = 1 


<1 = 1 

oo 

oo 


oo 

£ 

£ 

^' 2 ,2X” + 1 - 

• £ 

/*2=1 . 

h= 1 

* 

<2=1 

oo 

00 

* 

oo 

E %.*K 

£ 


• £ 

<*=> 

'*=i 


4=1 


£ z h.iK +k 


£ z i2 aK +k 


£ z i, k K +k 


(3.15) 


We are allowed to write (3.15) since D( 1, 1, . . 1), being a determinant, is the sum of a finite 

number (£!) of products of Ujj, and its asymptotic expansion as n^<x> is the sum of the products of 
the asymptotic expansions of the respective Ujj. By the multilinearity property of determinants, 
(3.15) is equivalent to 


H n~ £ • £ ( n z ip , p \ ( n k) ki. h \), 

/i = i <*— i \ p = i / \/?=i / 


(3.16) 


where K(£ 0 , lb . • I*) is the Vandermonde determinant defined by 


^(£o» • - * t f A t) = 


1 1 

lo li 

t 2 * 2 
€ i 


1 

I A: 


^0^1 ^A: 


(3.17) 


.Jw 


Since K(| 0 , . . £ k ) is odd under an interchange of the indices 0, 1, 

the appendix to this work, (3.16) can be expressed as 


. k, by Lemma A.l given in 





*/i.i 

Zi 2 ,l ■ 

’ • Z i k , 1 





OO 

H n ~ £ 

1 < /j < / 2 < . . 

•<<* 

*i.,2 

Zi 2,2 • - 

■ • «/*, 2 



(3.18) 




z i\,k 

z k,k • ■ 

• Z i k ,k 
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By (3.12), the most dominant term in the summation on the right side of (3.18) as w — oo is that for 
which /*! = 1, / 2 = 2, . . /*=£, provided that the determinant 


*1,1 

Z2,\ ■ ■ • 

Zk, 1 

* 1,2 

Z2,2 • ■ • 

Zk, 2 

Z\,k 

z 2,k • • • 

Zk,k 


(3.19) 


is nonzero. But since z,= (X/— l)u/, />1, we have 



(3.20) 


where Fis as defined in (3.11). Since /VO by assumption, FV 0 too; hence the first part of the 
theorem follows, with 


£>(1,. . .,1)= |^n (Xf-l^F^nX^Kd,^,. . .,X^)[l+o(l)] as n — «». (3.21) 

For the proof of the second part we proceed similarly. By (2.11), (3.3), and (3.10), we have 


1 00 


00 



00 



! E 

®/o X l 0 

E 

v u 

, x "„ +1 ••• 

E 

V <0 

■k:* 

! 4=i 


4 = 1 



4 = 1 



00 


oo 



oo 



E 

*/i.i x 7, 

E 

z ‘ i 

,1^ +1 --- 

E 

*<•. 


4 = 1 


4 = 1 



4 = 1 



00 


00 



00 



E 

Z '2,2\” 2 

E 

z h 


E 

Z/2 


4 = 1 

• 

<2=1 

• 


4 = 1 

• 


00 

• 

OO 

■ 


oo 



E 

z ‘k,k^” k 

E 

z ‘k 


E 

Z h. 


4=1 

4 = 1 


4=1 




Again from the multilinearity property of determinants 


J n~ E • • • E v io { IT z i p ,p \ ( FI • • •> 

<o=i /* = l \/>= l / \p = 0 V 


(3.22) 


(3.23) 


By Lemma A.l given in the appendix to this work, (3.23) can be expressed as 



(3.24) 




v h 

v h • • • 

v k 




*/i.i • • • 

z kA 

00 





G n~ £ 



z i\j 2 • * * 

Z ‘k |2 

1 ^ /"o < i\ . . 


. 

. 

. 



z io,k 

z h>k • • ■ 

Z ‘k,K 


\'o’ ’ * *> 


(A 1 *)- 


By the assumptions made following (2.11), there is a finite number of X; with moduli equal to |X*+ 1|. 
Let jX* +1 |=. . . =|X^ + ^>|X^ +r+ i|. From this and (3.12), it follows that the dominant term on the 
right side of (3.24) is the sum of those terms with indices /q= 1, /*i =2, . . ., i k -\ =k, /*=£ + £, £= 1, 2, 
. . ., r, that is, 

D(x n —s, . . ., x n + k — s) = ^ ^ ^£ + ^(Xj, • • *> ^k> ^k+t) 

\P = l / e=l 


V\ • 

• • V k 

v k+l 

Z\ y \ . 

• • Zk, 1 

Zk+e,i 

*1,2 • 

• • Z k ,2 

Zk+ e,2 

Z\ ,k • ■ 

■ Z k ,k 

Zk+l,k 


[1 + 0(1)] as oo. 


(3.25) 


Now the cofactor of in the determinant in (3.25) is F t which is nonzero since F^ 0. Therefore, 
for n suficiently large, the coefficients of t^+i, . . ., Vk+ r are nonzero. Since we also assumed that 
the Vj are linearly independent, the summation in (3.25) is never zero. Combining (3.21) and (3.25) in 
(3.7) results in (3.13). If |X^ + \\>^k+i\> then r= 1. In this case (3.14) follows from (3.21), (3.25), and 
the fact that 


^(io* •••»£*) = IT (£/-£/)• 

0 <i<j<k 


(3.26) 


Note: The condition on the determinant given in (3. 1 1) has some interesting implications when the 
normed vector space B is a complete inner product space. In this case, for each Qj in the dual space of 
B, there exists a unique vector qj in B, such that Qj(z) = {Qj,z) for every z in B> where (*,*) is the 
inner product associated with B. Then (3.11) becomes 


F= 


(QuV\) (Qi>v 2 ) . . . (quv k ) 
(42>V\) (? 2 ^ 2 ) * • • (Ql>V k ) 

( QbVx ) (q k ,v 2 ) • • • (Qk>Vk) 


* 0 . 


(3.27) 
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One of the consequences of (3.27) is that both sets of vectors Q k = Qk\ and V k — {i/j, . . 

Vfc) have to be linearly independent. Another consequence is that the intersection of the subspace 
span Q k with that orthogonal to span V k must be {0}. 

The asymptotic error analysis of the MMPE as given in Theorem 3.2 leads one to the following 
important conclusions: 

(1) Under the conditions stated in the theorem, the MMPE is a bona fide vector accelerator in the 
sense that 


Ikuzjjl = n 1” f bl£±l \ n 

Xn+k+1 -j|| |_\ h ) 


as/z-> oo. 


(3.28) 


This means that if x n —s, as /z — oo, that is, |Xj|< 1 , then s n> k~*s as /z — oo, and more quickly. Also if 
lim x m does not exist, that is, (Xj|>1, then s n ^s as /z — oo, provided that |A*+i|< 1. The reason that 

oo 

we write x n+ k+ \ in (3.28) is that s n> k in the MMPE makes use of the k + 2 vectors x n , x n + \, . . ., 

x n+k+ 1* 

(2) The result in (3.13) shows that when the MMPE is applied to a vector sequence generated by 
using the matrix iterative method described in Section 2.2, with the notation therein, it will be 
especially effective when the iteration matrix A has a small number of large eigenvalues (A: -many 
when s nt k is being used) that are well separated from the small eigenvalues. 

(3) By inspection of T(/z) in (3.13) and (3.14), it follows that a loss of accuracy will take place in 
s n,k when Xj, . . ., A* are close to 1, since ||T(/z) || becomes large in this case. When the vector 
sequence is obtained by solving the linear system of equations given in (2.12) by the iterative 
technique in (2.13), this means that if A has large eigenvalues near 1, there will be a loss of accuracy 
in s n k- In fact eigenvalues near 1 would cause the system in (2.12) to be nearly singular. 


4. Stability of MMPE 

Let us denote the yj of the previous section by y\ n>k ^ . Then the propagation of errors introduced in 

k 

the x m will be controlled, to some extent, by L (*•*> |; the larger this quantity, the worse the error 

y= o J 

propagation is expected to be. With this in mind, we say that s Hy k is asymptotically stable if 


k 

sup \yj n ’ k) \<oo. 
n j = o 


(4.1) 


Since Jj = 1 by (3.2), then Lto*) | >1, so that the most ideal situation is that in which 
j-o j = o 

y(n,k) >o, 0<j<:k, for n sufficiently large. The following theorem shows that for the type of 
sequences considered in Theorem 3.2, s n> k as obtained from MMPE is asymptotically stable, and that 
y(n,k) > 0 ^ o<y<£, for sufficiently large n, whenever X„ 1 </</:, are real and negative. 

Theorem 4.1: Under the conditions stated in Theorem 3.2, s Ht k is asymptotically stable. 

Proof: By (3.5), it is sufficient to show that yj n,k \ 0 <j<k, stay bounded for /z — oo, which in turn 
guarantees (4.1). Now 



«/U • 

• u n+j- 1,1 

u n+j+ 1,1 • 

■ u n + k,l 

Un.2 • 

• Un +7-1,2 

M n +j + 1 ,2 • 

■ u n + k,2 

u n,k * • 

• u n+j-\,k 

u n+j+l,k • • 

• u n + k,k 


(4.2) 


Substituting the asymptotic expansions of u m j as given in (3.10) and using the multilinearity property 
of determinants, we obtain 


OO 

Nj~ L ■ 

ii = \ 

oo / k \ / k \ 

• • e ( n ztp.pj ( n x ? 

ik ~ 1 \P=1 / \P=1 / 

(cA.,. 

. Xf t ) &sn-<x>. 

(4.3) 

where 


i a 1 

A 1 •• 






1 «2 ••• A' 

^2 +1 •• 

• 



wax. ■ ■ ■ 

. **>=(- o' 

A' 

A' '• 

• ^ 


(4.4) 


Since Cj(%\ 9 . . £*) is odd under an interchange of the indices 1, . . M k. Lemma A.l in the 

appendix again applies, and we have 




^/|,1 1 

* * * ^’*,1 



OO 

Nj~ L 

1 </*! </' 2 <. . 

.<i k 

^/|,2 £/ 2 »2 

• ■ • ^,2 

Cj(\, • • 

-K)(nK 


1 


* • • 



By (3.11), (3.19), (3.20), and (3.12) 



Nj=FCj(\ u . 


'(£.*) 

1 [l+o(D] 

as n — oo. 



(4.5) 


(4.6) 



Combining (4.6) and (3.21) in (3.5), and using (3.20), we obtain 


Obviously, (4.7) also implies that [ 7 >J n - k ) |<oo for sufficiently large n. This then proves (4.1). 
Inspection of (4.4) reveals that C\(X 1 , . . X*) is the cofactor of X 7 in the first row of 


1 X . . . X* 


1 X. . . . X* 


K(X, X,, • . ., X*) = I. 


1 X*. . . . K 


that is, 


K(X, X„ . . ., X*) = £ C y (X b . . ., \ k )\l 
7=0 

Combining (4.7) and (4.9), we obtain the following interesting result: 


fa < 410 > 

n-ooy = 0 7 Kl.X,,. . ..X*) 

Invoking (3.26), (4.10) finally becomes 

lim t rj’ M V= n (t=t)- (4.11) 

n-+ OO y = 0 1=1 ' A/ 1 / 

Since K(X, Xj, . . ., X^) is a polynominal of degree k and X and vanishes when X = X/, 1 </<£, we 
have 


K(X, X lt . . ., X*) = C*(X,, . . ., X*) A (X— X p ). 

/7=1 


(4.12) 


Upon expanding the product on the right side of (4.12), and comparing with (4.9), we get 


C*-y(Xi, ■ . x*) =(- lyC^Xj, . . ., X*) £ \ ■ ■ ■ H> 1 -j- k - (4.i3) 

1 <i { <. . ,<ij<k 
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It is obvious from (4.13) that if X/, 1</<A:, are real and negative, then sgn C^_y(X 
X*) =sgn C k (\ u . . X*) for 0 </<£, and this implies that >0, 0<y <k, for n sufficiently 

large. If X/, 1 </<£, are real and positive, then (4.13) implies that yj n > k )yj"*p <0, 0</‘ <&- 1. 

5. Convergence and Stability of TEA 

In this section we shall consider the convergence and stability properties of the TEA, which is 
obtained from Approach (a) or (b) in conjunction with Method (3) of Section 2.2. The TEA can be 
summarized as follows: The approximation s n k to s is given by 

k 

s n t k =z 2 yi x n + i» (5-1) 

i—0 

where y / are obtained from the equations 
k 

E y> = 1 

i = 0 

k 

2 j 7/G(«m+i) =0, n<m<n + k-\ y 

i*-0 

with w m = Ax m =x m+ j —x m , m> 0, as before. When 7^^0, equations (5.2) are equivalent to (2.18) 
and (2.22), in which dj=Cj, w m+i = u m+i> 0</<&-l, and w m = —u m+k , as can be verified by 
inspection. 

We now write the equations in (5.2) in the form 
k 

t 7/= 1 
i = 0 

k 

E 7i u n + ij = 0, 

i = Q 

where this time u m j = Q(u m+ j_\), 1 </</:. Defining D(o$, a\ 9 . . ., o k ) as in (3.3) but with 

the u m j of Section 3 replaced by the new u m j y we see that s ntk for TEA is given exactly by (3.6), as 
can be verified by Cramer’s rule. Let us also define z/j=Q(z,)X(~ 1 = Q(u/) (X/- l)Xj“ /> 1, 
1 <y <k, where the z \ are as defined in Section 3. Then under the assumptions preceding (3.10), (3.10) 
holds with the Zjj of Section 3 replaced by the new zij . 

If one follows the proof of Theorem 3.2, one realizes that (3.18) and (3.24), which form the most 
important parts of it, are consequences of (2.11) and (3.10). Consequently (3.18) and (3.24) remain 
true for the TEA provided the Zjj of Section 3 are replaced by those of the present section. Starting 
with these observations, we now prove the following theorem: 

Theorem 5.7: Assume that 

Q{vd* 0, 1 </<*, (5.4) 

and that all the conditions of Theorem 3.2 (with the exception of (3.11)) are satisfied. Then, for all 
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k 

sufficiently large n, Nj^ 0, hence s n ^ as given in (3.6) exists. Furthermore, 
j~0 

^Ar— *s= A(/i) X" + j[l +o(i>] asn-oo, (5.5) 

where the vector A(w) is nonzero and bounded for all sufficiently large n. If, in addition, 

l\t+ll>l\t+2|. then 


Vl 

v 2 . 

• v k+l 





*1,1 

*2,1 • 

• Zk+ 1,1 



k 

n (\t+i - ^/) 


* 1,2 

2,2 • 

• Z/t+1,2 



i=i 


V(\i, . 

* •> ^k) 



• 


• 

*u 

Z2,k ■ ■ 

■ • Zk+\,k 

! 





Proof: The proof of this theorem proceeds along the same line as that of Theorem 3.2, using the 
additional relation 


*u 

*2,1 • 

■ Z k , 1 



*1,2 

*2,2 • 

• */t,2 

= K(X,.. . 

•.x*)[nG(«i)] 

*1 ,k 

*2 ,k • • 

• Zk,k 




□ 

Note: When .the normed vector space B is a Hilbert space, the condition (5.4) has the following 
implication: Let q be the unique vector in B for which Q(z) = ( q,z ) for every z in B. Then (5.4) 
implies that q cannot be orthogonal to any of the vectors Vj, 1 </<£. 

As a result of the asymptotic error analysis of the TEA given in Theorem 5.1, we can draw 
conclusions that are identical to those about the MMPE given at the end of Section 3. We find that 
(3.27) is replaced by 


lfou-g|| =n 



as /? — oo. 


(5.8) 


since s n k for TEA is formed by taking into account the 2A: -h 1 vectors x n , x n + . . ., x n + 2k> instead 

of the £ + 2 vectors x n , x n+ . . ., x n +k+\ used to form s nj k for the MMPE. This in turn implies that 
the MMPE is a more economical vector accelerator than the TEA, since it attains the same rate of 
acceleration as the TEA while using approximately half the number of vectors. 

Finally the stability properties of the TEA are very similar to those of the MMPE as is stated in the 
following theorem. 



Theorem 5.2: Under the conditions stated in Theorem 5.1 s Htk is asymptotically stable. 
Furthermore, (4.7) and (4.11) hold too. 

Proof: Similar to that of Theorem 4.1. 

□ 

The conclusions that were drawn from the stability analysis of the MMPE in Section 4, are true for 
the TEA too, as some analysis reveals. 


6. Numerical Examples 

In Section 3 we analyzed the convergence properties of s n>k for the MMPE as n — oo, and derived an 
asymptotic error estimate for it, obtaining at the same time its rate of convergence. In this section we 
apply the MMPE and the MPE to three vector sequences obtained as iterative approximations to 
linear systems of equations. The numerical results verify the conclusions of the asymptotic error 
analysis of Section 3. They also indicate that the MMPE and the MPE have very similar 
performances. 

In all the examples below the MPE is implemented by solving the generally overdetermined system 


k-1 

^ c i u n + i = ~ u n + k 
/ = 0 

/ k 

^2 c j> 0 </<£, with c k = 1, in 
J = 0 

(2.17). As before Uj=Xj+i-Xj, j = 0, 1, . . .. The MMPE, on the other hand, is implemented by 
solving the linear system of k equations 

k - 1 

Ti c » u n+iJ = ~Un + kJ> 1 ^j^k, 

7 = 0 

where u m j denotes the j th component of the vector u m . That is to say, we pick the linear functional 

Qj in (2.17) to be the projection operator onto the subspace spanned by the/h un it vector, j= 1 

k. 

/ k 

^ Cj, 0</<k, with c k =\, in (2.17). 

J - 0 

Example 1. The vectors x, are obtained by setting jc 0 = 0 and x i+ j -Axj+ b, where A is the iteration 
matrix associated with the Gauss-Seidel method for the system of linear equations Cx=d, where 


C= 


2 

1 

3 

.4 


1 

-3 

1 

5 


3 4 

1 5 

6 -2 
-2 - 1 _ 


and d and/or b are determined by requiring that the solution s to Cx=d be the vector with all its 
entries equal to 1. The eigenvalues of A are approximately X( = X 2 = - 2.3500 ±/2. 0506, 
X 3 = -0.0228, and X 4 = 0 . Therefore, lim Xj does not exist. In Table I we give the errors H^/t - s II® 

/— 00 

computed in the norm for k=2 and 0<ns5, both for the MMPE and the MPE. 



TABLE I. — 4. NORMS OF THE 
ERRORS s n> 2 —s FOR 
EXAMPLE 1, COMPUTED 
USING MMPE AND MPE, FOR 
k=2 AND 0</z<5 

[The numbers have been rounded 
to one significant decimal 
digit. The base iterations x j 
diverge.] 


n 

MMPE 

MPE 


Ikn,2-S|L 

lk/r.2 -s 1. 

0 

6x 10' 1 

1x10“ 

1 

8x10-3 

7x10-3 

2 

2x10-* 

2x 10 -4 

3 

4xl0-« 

4x 10 -6 

4 

1 x 10“ 7 

9xl0- 8 

5 

9xlO-io 

9x10-'° 


Example 2. The vectors x / are obtained by setting Xq= 0 and x l + 1 =AXj+b, i= 0, 1, . . where 


5 2 
2 6 
1 3 
1 1 
1 


A =0.06x 


1 1 
3 1 1 

6 3 11 
3 6 3 1 
13 6 3 
113 6 
1 1 3 

1 1 
1 


1 

1 1 

3 1 1 

6 3 11 

3 6 3 1 1 

1 3 6 3 1 

113 6 2 

112 5 . 


and b is determined by requiring that the solution s of the system x=Ax+b be the vector with all its 
entries equal to 1. The eigenvalues of A are all real and in (0,1), and are approximately Xi =0.8965, 
\2 = 0.7318, \3 = 0.5297, X 4 = 0.3600, . . ., Xu =0.0313, in decreasing order. Since Xj < 1, the 
sequence x it i=0, 1, 2, . . ., converges. 

In Figures 1 and 2 we give the results of the computations for ||s n> *— s H* using both the MMPE 
and the MPE with k= 1 and k= 2, respectively. The figures also include ||x n+ *+i -s ||oo> 




10 I— 


9 H 


8 r~ 


7 1 - 


2 4 
S’ 


8 


2h 


ih 


o 


4 


A -l°9io II x n+k+l " 

□ -log 10 II S n k -SlL for MMPE 

O -iog 10 I! S P( k - S| ro for MPE 

O 

O 

O □ 

o □ 
o □ 
o □ 
o □ 
o □ 
o □ 
o □ 



n + k + 1 


Figure 1. - Results for Example 2 taking k= 1. 
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A "lotJioll x n+k+l " Sllco 
[— □ -log 10 II S n k -S||„ forMMPE 

O 'log 10 II S n _ k - SIL for MPE 


8h 


S 5 


O 


I a I 8 1 A 


□ 


□ 


□ 


J 1 I L 


o 


□ 


□ 


J I L 


i 


10 12 14 16 18 20 22 24 26 28 30 

n + k + 1 


Figure 2. - Results for Example 2 taking k-2. 


Example 3. The vectors x t are obtained by setting x 0 =0 and x i+ 1 =AXj+b, i = 0, 1, . . where 


12 11 
11 11 10 
10 10 10 9 


A =0.04 X ’ ’ | ' 


2 2 2 2 ... 2 1 

1 1 1 1 ... 1 1 


and b is determined by requiring that the solution 5 of the system x=Ax+b be the vector with all its 
entries equal to 1. The eigenvalues of A are all real and are approximately Xi = 1.2892, X 2 = 0.8080, 
X 3 = 0.4924, X 4 = 0.2785, . . ., X 12 =0.0012, in decreasing order. Since Xi>l, the sequence Xj, i= 0,1, 
2, . . ., diverges. In figures 3 and 4 we give the results of the computations for s ||oo using both 
the MMPE and the MPE, with k=2 and k= 3, respectively. 



- Iog 10 Uoo norms of errors) 


10 r~ 


9 — □ - log 10 i| S n> k - SL for MMPE 

° “log 10 II S n> k - SlL for MPE 

8 — 


7 — 


6 


5 — 
4 — 

3 — 
2 — 



O 

□ 


1 — 

0 — 
-1 

0 


§ 

0 

□ 

□ □ 

o o o 

1 I I I I I I I I ! l_l J I I 

2 4 6 8 10 1 2 14 16 18 20 22 24 26 28 30 

n + k + 1 


Figure 3. - Results for Example 3 taking k = 2. The base iterations Xj diverge. 


21 


on 



10 


□ -log 10 l S n>k -SIL forMMPE 
O -log 10 ll S nk - Sll* for MPE 


O 


□ 


E 5 


o 


□ 


1 \— 


0 I— 


O 

□ 


□ 

o 


0 2 4 6 


I I I l l i I l i i I I 

8 10 12 14 16 18 20 22 24 26 28 30 

n + k + 1 


Figure 4. - Results for Example 3 taking k — 3. The base iterations Xj diverge. 
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Appendix 

Lemma A. 1: Let / 0 , i u . . i k be integers greater than or equal to 1 , and assume that the scalars v ioi 

i k are odd under an interchange of any two indices / 0 , ...» /*. Let 07 , i > 1, be scalars (or vectors), 
and let uj, i* i, i ^j<k, be scalars. Define 


and 


Jfc,N= 


TV . 

• E <% 

<*=i > 

< A: 

n '< P 

<p=i 

■pj v ‘o. ■ ■ 




% 

a h . . . 




'<0.1 

'<i,l • • • 

'<*.1 

E 


'<0.2 

'<1.2 • • • 

'<*.2 


1 <• • • <i/c^N 


,k U u k • • * U k ,k 


(A.l) 


> • • •> 


(A. 2) 


where the determinant in (A. 2) is to be interpreted in the same way as D(o o, . . o k ) in (3.3). Then 

h,N~Jk,N • (A. 3) 

Proof: Let L k be the set of all permutations of {0, 1, . . ., k}. Then by the definition of 
determinants 


Jk,N ~ £ £ ( s & n *>%()) ( IT * • - *k' 

l<lo</‘i<. . .</*<N tCE* \p=l / 


(A. 4) 


Now 7 i ^ 7 r (/?) == /?> and sgn tt " 1 =sgn ir for any permutation 7 r€E*. Hence 

J k,N= E E ( s 8 n 1r-1 )% 0) ( n 'x-'Tfp)' 

l£/ 0 </,<. . ,<i k ^N ir€E* \ /7 = 1 / * W T 


(A. 5) 


By the oddness of v iu m m m / , we have 

'Tr-^CO)’ ’ * -’'n-Mp) 


(sgn 7T- 1)^/ , / ... =Vi 

— — 1 T 


(A. 6 ) 
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Substituting (A. 6) in (A. 5), we obtain 


J k,N- £ 

1 </ 0 </ I <. . . </*< N 



(A. 7) 


Since Vj Qi / is odd under interchange of the indices /q* . . ., it vanishes when any two of these 
indices are equal. Using this fact in (A. 1), we see that I/c 3 n is just the sum over all permutations of the 
distinct indices /q, . . i k . The result now follows by comparison with (A. 7). 

□ 

Note that when the a,- are scalars, (A. 3) remains true also for the case in which the U/ 0 , . . i k are 
vectors. When the and v- lQi # m i k are vectors, (A. 3) still holds provided a/ 0 V/ 0t . . i k is interpreted as 
a direct (tensor) product. 
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